% Calculate mean WTP to avoid Newberry site closure by Newberry hunters
% only

clear

rng(1196)

% Load data, assign parameters
load('sigma1save_myopic.mat')
load('sigma0save_myopic.mat','phi0', 'Vb0', 'zeta0');
% pit = [thetastar(4) thetastar(5) 1-thetastar(4)-thetastar(5)];            % Type probabilities
s = thetastar(end);                                                         % Cognitive cost parameter
mu = thetastar(1);                                                          % Travel cost parameter
eta = [thetastar(2) thetastar(3) 0];                                        % Type fixed effects
chi = [0 thetastar(6:end-2) 0];                                             % Hunt baseline utilities
phi0(end,:) = zeros(1,K);
phi1(end,:) = zeros(1,K);
tc = [zeros(size(tc,1),1) tc zeros(size(tc,1),1)];

% Specify type, travel cost distributions from model estimates and data
typedist = makedist('Multinomial',pit);                                         
tcdist = makedist('Multinomial',PTC);

% Initialize result vectors, set number of iterations
iter = 1e5;
tctracker = zeros(iter,1);
typetracker = zeros(iter,1);
pptracker = zeros(iter,1);
jstar0tracker = zeros(iter,1);
jstar1tracker = zeros(iter,1);
phiv0 = zeros(iter,1);
phiv1 = zeros(iter,1);

% Run simulation
parfor idx1 = 1:iter
    
    % Draw travel cost profile, type
    tcdraw = random(tcdist);
    tctracker(idx1) = tcdraw;
    typedraw = random(typedist);
    typetracker(idx1) = typedraw;
    
    % Define pref pt distribution conditional on type, tc profile, then
    % draw pref pt stock
    ppdist = makedist('Multinomial',zeta0(:,:,tcdraw,typedraw)');
    ppdraw = random(ppdist);
    pptracker(idx1) = ppdraw;
    
    % Define type-1 EV and exponential error terms (from hunting (e) and cognitive cost (w))
    e = -log(-log(rand(H-1,1)));
    w = log(1 - rand())/-s;

    error0 = phi0(:,ppdraw).*[0; e; 0] + [zeros(H,1); w];
    error1 = phi1(:,ppdraw).*[0; e; 0] + [zeros(H,1); w];
    [~,jstar0] = max(Vb0(:,ppdraw,tcdraw,typedraw) + error0);
    [~,jstar1] = max(Vb1(:,ppdraw,tcdraw,typedraw) + error1);
    jstar0tracker(idx1) = jstar0;
    jstar1tracker(idx1) = jstar1;
    
    % Calculate flow utility before, after closure
    phiv0(idx1) = phi0(jstar0,ppdraw)*(chi(jstar0) + eta(typedraw) ...
        + mu*tc(tcdraw,jstar0)) + error0(jstar0);
    phiv1(idx1) = phi1(jstar1,ppdraw)*(chi(jstar1) + eta(typedraw) ...
        + mu*tc(tcdraw,jstar1)) + error1(jstar1);
        
end

WTP_Newberry = mean(phiv1(jstar0tracker==17|jstar0tracker==18|jstar0tracker==19) ...
    - phiv0(jstar0tracker==17|jstar0tracker==18|jstar0tracker==19))/abs(mu/200)

WTP_Not_Newberry = mean(phiv1(jstar0tracker~=17&jstar0tracker~=18&jstar0tracker~=19) ...
    - phiv0(jstar0tracker~=17&jstar0tracker~=18&jstar0tracker~=19))/abs(mu/200)

Newberry_prop = sum(jstar0tracker==17|jstar0tracker==18|jstar0tracker==19)...
    /size(phiv1,1)

% save('Results.mat')
